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Abstract 

We present a short overview over recent MD simulations of systems of fully 
flexible polyelectrolyte chains with explicitly treated counter ions using the full 
Coulomb potential. The main emphasis is given on the conformational properties 
of the polymers, with a short discussion on counter ion condensation. 

1 INTRODUCTION 

Polyelectrolytes represent a broad and very interesting class of materials [|TJ that enjoys an 
increasing attention in the scientific community. Even though the theory of neutral poly- 
mer systems is well developed, polyelectrolytes remain one of the least understood states 
of the soft matter area. Theoretically this is mainly due to the long range Coulomb inter- 
action. Simple scaling theories, which have been proven so successfully in neutral polymer 
theory, have to deal with the additional length scales set by the electrostatic interaction. 
It is the delicate interplay between electrostatic interaction, and the conformational de- 
grees of freedom, which in turn are governed by a host of short range interactions, which 
proves to be a difficult question for the theory. This results in essence in the fact that 
the analytical theory is only able to treat two limiting cases, namely the case of high salt 
excess, resulting in effectively screening down the electrostatic interaction to treat it as 
a perturbation, or the case of an overwhelming dominance of the Coulomb force, which 
results in a strongly elongated chain. Unfortunately it is just the intermediate case, which 
proves to be the most interesting regime in terms of application, experiment and theory. 

How difficult it is for the analytical theory to treat this intermediate regime could be 
seen in the recent controversy over the correct exponent y, which governs the variation 
of the electrostatic persistence length L e with Debye screening length r^.. The theory is 
already simplified because instead of the full Coulomb interaction it assumes a screened 
Debye-Hiickel potential Uonir) — AfffcBT exp( - r//rg ' > , where A# is the Bjerrum length given 

2 

by Xb = 47re \ ksT i with e denoting the elementary charge and e s ,e r are the dielectric 
constants. The two competing theories both predicts L E oc rf>, however the exponents 
are y = 1 or y = 2, depending which perturbation approach one believes more. To test 
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the theoretical approaches, a series of molecular dynamics (MD) and Monte Carlo (MC) 
simulations have been presented with the result, that neither theoretical approach is 
correct in the intermediate regime, but that the exponent y can vary continuously between 
0.5 and 2. Only very recently a field theoretic renormalization group analysis was able to 
shed some additional light onto this problem||. 

A serious drawback of the Debye-Htickel theory is that it presupposes only weak vari- 
ations in the density of the counterions, and neglects all correlation effects. One can 
therefor doubt if this potential is able to describe the chain conformations adequately 
Q. Also counterion condensation can not properly be handled. To tackle these questions 
computer simulations provide the necessary tools, because here one does not need to rely 
on excessive simplifications. One can treat the counterions explicitly and can also utilize 
the full long range Coulomb potential. We try to summarize here recent results of MD 
simulations on systems of polyelectrolyte chains in solution. Details of the studies can be 
found in the original literature [H, []. 0] 



2 MODEL 

Our model of a flexible polyelectrolyte chain consists of N p bead-chain polymers with 
N m monomers which are located in a simulation box with periodic boundary conditions 
(3D torus). From these monomers a fraction / is monovalently charged. The number of 
counterions, N c , with valence v is then chosen such that the overall system is electrically 
neutral. If we have, for example, N q charges on a polymer chain, then the number of 
counterions with valence v is given by N c = NpNq , and the total number of charges in the 
system is N tq = 2N C . All hydrophilic monomers are given an effective size through a pure 
repulsive Lennard- Jones potential, representing thus a polymer chain in a good solvent: 
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All chain monomers, which model the hydrophobic or poor solvent case, interact via a 
standard Lennard- Jones potential with attractive part: 
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All chain monomers are in addition connected by the FENE (finitely extended nonlinear 
elastic) bond potential, 
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Unsssir) = --kB*ln\l-j^). (3) 
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All charged monomers posses in addition the full Coulomb energy 



EM,) = l \ (4) 



where v i is the valence of the i th charged monomer. 

The electrostatic energy of the box is calculated with the Ewald formula 

E c = E ir) + E {k) + E (s) (5) 

where the contribution from real space E^ , the contribution from reciprocal space 
and the self energy are given by 

(r) = ly y eMa\r ij + mL\) 

2 ^ ^ H H3 r« + ml v ; 

E(k) = ^E^ e " fc2/4Q2 ip( k )i 2 ( 7 ) 

= -^=5> 2 (8) 



and the Fourier transformed charge density p(k) is defined as p(k) = Ylj=i Ij e ~ % k Fj • 

The Ewald parameter a, which has the dimension of an inverse length, tunes the 
relative weight of the real space and the reciprocal space contribution, but the final result 
is of course independent of a. The k- vectors form the discrete set {27rn/L : n 6 Z 3 }, and 
we use tin-foil boundary conditions, hence there is no dipole correction term. The standard 
Ewald method has computational effort at best of 0(Nf ( j 2 ) J which limits the size of the 
samples. To overcome this limitation, we use for E^ a particle mesh Ewald algorithm 
(PME) that uses discrete charge assignments and fast Fourier transforms (FFT), and 
accordingly scales for E^ as 0(N tq (\og N tg ) 3 ^ 2 ) [§]. Note, however, that this is not the 
most effective algorithm, as was shown recently 0. 

In all our performed molecular dynamics (MD) simulations we had no salt ions and 
no explicit solvent molecules. The sol vent, however, was implicitly taken into account 
via a Langevin thermostat with damping constant T = r _1 , with time step 0.00125r at 
constant temperature fcgT = le. The number of MD steps was chosen such that the 

typical observables like the end-to-end distance Re = y (Re) or the radius of gyration 
Rg had sufficiently relaxed, which happened usually after 500 000 up to 5 000 000 MD 
steps. 

3 STRONGLY CHARGED POLYELECTROLYTES 
IN GOOD SOLVENT 
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3.1 MONOVALENT COUNTERIONS 



The first investigation of totally flexible many chain polyelectrolyte systems in good sol- 
vent with explicit monovalent counterions was performed already some years ago f| . The 
simulations were done mostly with systems of 8 or 16 chains with N m = 16, 32, and 
64. Instead of the PME algorithm a spherical approximation in a truncated octahedral 
simulation box was used, which is for smaller values than N tq ~ 500 faster than the PME 
method. More details of the whole study can be found in ||. 

In this work it was demonstrated that the simple bead-spring model can actually 
be compared to real polyelectrolytes, because the experimental values for the osmotic 
pressure and the maximum position in the interchain structure factor are successfully 
being reproduced. One of the important findings was that the rodlike chain is found to 
be a rarity already in the dilute limit. Counterion condensation can dramatically shrink 
the polyelectrolyte chain. The end-to-end distance shortens significantly as the density 
increases from the dilute saturation value to the overlap value. The chain structure is 
highly asymmetric at the very dilute saturation density and the scaling with respect 
to N m is asymmetric, but as the overlap density is approached, the structure is less 
asymmetric and the scaling becomes approximately symmetric. On long length scales the 
chain structure continuously changes from very elongated to neutrallike coils. Yet, on 
short length scales, the chain structure is density independent and elongated more than 
neutral chains. The findings for the single Bjerrum length Xb = 0.883a are summarized in 
Fig. p]. It was found that in the dilute limit the scaling for the extension perpendicular to 
the chain was R± oc jV a65-a7 °, and for the extension parallel it was R» cx jV a90_L00 . Near 
the density, where the rodlike chains in disordered solution would overlap, p ~ iV -2 , R± 
starts to grow on the expense of R» until at the overlap density p* the effective exponent 
is about 0.82. This transition regime ranges from p ~ N~ 2 to about p ~ iV~ 1,4 where the 
coils start to overlap and one eventually reaches v — 1/2 in the semidilute regime. The 
exponents reported should not necessarily be taken as asymptotic (N — > oc), however 
they should be relevant for many experimentally systems. 

In a recent study a large MD simulation of 200 chains with N m = 32 was performed 
to determine if with these large systems (N tq = 12 800) any large scale ordering phenomena 
could be observed. The idea of the onset of ordering is one of the possible explanations for 
the Fuoss-Strauss effect in the reduced viscosity rj r , which shows a dramatic increase and 
subsequently a decrease of r] r , as the polymer concentration is decreased. The pure MD 
simulation showed so far no sign of correlations or ordering, however it was found that 
the center of mass diffusion of the monomers was just too slow to move the monomers far 
enough. It needs definitely accelerated algorithms (MC or MD) to settle this question in 
a large scale computer simulation, which will be pursued in the future. 
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Figure 1: A structure diagram for A# = 0.833cr, which shows the dilute limit, the con- 
traction below the overlap density p*, and the completely screened limit at large p. 



In this set of first simulations we usually considered eight chains in the central box with 
a chain length of N m = 106 and a charge fraction of / = 1/3. We varied the valence 
of the counterions from one to three, and chose always Xb = 3a, thus the monovalent 
systems were just below the critical Manning threshold. The weighted density of charged 
particles, p c = ^ * ? was varied from 10 -1 cr -3 . . . 10~ 6 o~~ 3 . The starting configurations 
were set up as rod configurations with randomly distributed counterions, or as a random 
walk, if the rods would not fit in the central box. That both initial configurations lead to 
the same equilibrium value for Re can be seen in Fig. |^. Another point to note is, that 
both time series converge very rapidly, even before they are fully equilibrated. This is 
due to the fact, that the Coulombic repulsion of the chain monomers provide a very fast 
initial rearrangement of the chain, but that the equilibrium conformation is only reached, 
after the counterion cloud has itself arranged around the polymers. If one looks on the 
end-to-end distance Re as a function of density for the different valences, see left side of 
Fig. H], one notices that for high polymer concentrations all three valences have basically 
the same Re, they show SAW behavior. For the semi-dilute range the monovalent systems 
elongate very fast, whereas the divalent systems elongate much slower, and the trivalent 
system actually contracts. This is due to the fact that a trivalent counterion can bind more 
than one charged chain monomer, so that the chain can "wrap around" the counterion, 



3.2 MULTIVALENT COUNTERIONS 
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Figure 2: Time series of the end-to-end distance Re versus simulation time t/r, plus a 
typical equilibrium configuration for a trivalent system at p = 10 _5 o~~ 3 . 



or the counterion can even bridge two far apart regions of the chain, as can be seen in the 
configuration in the upper part of Fig. For lower density the trivalent chains elongate 
much slower, so that they will reach a rod like configuration at much lower concentrations. 
Even at this dilution the chains still show much structure due to the discrete counterions. 

Another interesting quantity to look at is the fraction of counterions found up to a 
distance r, which we denote by the integrated ion distribution function P(r), shown on 
the right side of Fig. |3|. Especially for low concentrations, this function gives a good 
estimate to determine the fraction of condensed counterions, because one can clearly see 
a plateau value for P(r), meaning that the ions can be separated in bulk counterions 
far apart from the polymer, and a distinct layer of condensed counterions. For higher 
concentrations and lower valences, however, this distinctions gets less pronounced. One 
observes that the condensed fraction of counterions, which we denote by C, depends 
strongly on the density and valence. With increasing density and valence the condensed 
fraction increases, and assumes for the highest investigated density p = 10 _1 <7~ 3 the value 
C = 1 for all three valences. The Manning predictions for these systems are C = 0(v = 1), 
because for the monovalent systems we are slightly below the critical Manning parameter, 
and C ~ l/2(t> = 2) and C ~ l/3(u = 3) for the other valences. 

A third interesting aspect to look at is the counterion exchange dynamics. If one 
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Figure 3: Left: Re vs. p for valences v = 1...3; right: integrated ion distribution function 
P(r) vs. distance r for p = 10~ 3 a~ 3 and p = 10~ 6 cr~ 3 , each for valences v = 1...3. 

defines the set of condensed ions as IK^(t) , which contains all ions being within distance d 
on polymer j at time t, then one can define a conditional condensation fraction Cd{{t\to) 
as follows: 

\ K(i„)| l j 

This equation simply means, count the number of ions which are on time t as well as on 
time to condensed on polymer j, normalize by the number of condensed ions on polymer 
j at time to, and average then over all polymer chains j. The normalization is such, that 
c,i(^o|^o) = 1- Fitting this quantity to an exponential Ansatz Cd(t\to) ~ A + B exp — t/t re iax 
gives information about the exchange dynamics of the counterions. We find that with 
increasing valence at constant density the mobility decreases, and that the mobility also 
decreases with decreasing density at constant valence. The value of t re iax is of the order 
of 10 OOOr, showing the very slow dynamics at low densities. 

4 POLYELECTROLYTES IN POOR SOLVENTS 

The investigations so far dealt only with polyelectrolytes in good solvent, to facilitate 
the comparison to known systems such as neutral polymer systems. Most experimentally 
investigated polyelectrolytes, however, posses a hydrophobic backbone. This results in 
a competition between the solvent quality, the Coulombic repulsion, and the entropic 
degrees of freedom, which can lead to totally different behavior 0. The simulations pre- 
sented below used 16 chains of length N m = 94, with a charge fraction of / = 1/3, and 
monovalent counterions. The hydrophobic interaction strength was tuned by means of 



the Lennard- Jones parameter elj of Eq. (0) that was chosen such that the finite chain 
shows an effective random walk behavior (elj = 1.5) for the density p = 10~ 4 <7 -3 . 

The polymer density p can be used as a very simple parameter to separate different 
conformation regimes. This can already be seen in the plots of the end-to-end distance 



R P and r 



versus p in Fig. f|. At very high densities the electrostatic interaction 



is highly screened, so that the hydrophobic interaction wins, and the chains collapse to 
dense globules. If one slightly decreases the density, the chains can even contract further, 
because there are no more steric hinderences from the other chains or counterions, and the 
screening is smaller. The collapsed globules, however, have still a net charge, and repel 
each other, so that this phase resembles a charged stabilized colloid or microgel phase. 
With decreasing density the electrostatic interaction will dominate over the hydrophobic 
one. The chains will tend to elongate, assuming pearl-necklace conformations, see Fig. [5|, 



as they have been predicted for weakly charged polyelectrolytes in Ref. []T0| . The more 
the chain stretches, the smaller become the locally compact regions. Reducing the density 
even further then leads to a strongly elongated blob pole, which is however, due to the 
intra chain entropy and the hydrophobicity, still far from an ideal stiff cylinder, and shows 
local structure. 
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Figure 4: Re (left) and r (right) versus density p for hydrophilic (phil), weak hydrophobic 
(clj = 0.5), and strongly hydrophobic (e^j = 1.5) chains 



5 CONCLUSION 

We summarized the recent efforts which were pursued to learn more about linear flexible 
polyelectrolytes in solution. The investigated systems are already large enough to repro- 
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Figure 5: Typical polyelectrolyte conformation for the densities p = 2 ■ 10 4 ct 3 (left) and 
p = 2 ■ lCr 5 (j~ 3 (right). Only counterions within 10 a of the chain are considered. 



duce experimental measurements such as the osmotic pressure or the intrachain structure 
factor. The incorporation of multivalent counter ions alters the chain structure drastically. 
Such aspects can clearly never be captured by a Debye-Hiickel approach. The integrated 
ion distribution function P(r) seems to be a good measure of how one can separate the 
condensed layer of counterions from the bulk counterions, and shows clearly that the 
fraction of condensed ions is density dependent. 

With the beginning investigations of polyelectrolytes in poor solvents one can also 
compare the simulations more easily to experiments. An interesting aspect of hydropho- 
bicity is, that already weak hydrophobicity can alter the chain structure, and that with 
strong hydrophobicity the system can produce a physical gel. The behavior of Re versus 
the density looks very similar to the behavior which is found for hydrophilic systems with 
trivalent counterions. 

Still many things remain to be done. To name just a few, the whole range of semi- 
flexible linear chains is still to be investigated. A detailed study of Manning condensation 
is also lacking. To analyze properly any sort of aggregation phenomena, one needs to sim- 
ulate much larger systems. This can, up to now, only be pursued, if one has an efficiently 
parallized code, which incorporates a mesh Ewald method or any other algorithm whose 
CPU-time scales linearly with the total number of charges N tq . 
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